function G=initialize_metric_QSE(xiIn)
%initialize_metric_QSE -- initialize structured metric in QSEPACK.
%   Given a cell array xiIn containing algebraic state coordinates,
%
%      G = initialize_metric_QSE(xiIn);
% 
%   initializes and returns a structure G specifying the metric matrix
%   in a compressed representation that is suitable as an argument for 
%   the QSEPACK routine compute_product_QSE(...).
%
%   The input xIN is a cell array of complex row vectors which have been
%   initialized according to the following indexing convention (specified 
%   in LaTeX, see the QSEPACK documentation for details):
%
%      xiIn{j,k}(1,l) = ${}\sp{j}\xi^{k}_{l}$
%
%   Note: for given k, all xiIn{:,k} vectors must have the same length.  
%   Violating this condition generates a run-time error message.
%
%   An alphabetized list of all the elements of the G structure follows:
%
%   G.coefs   = [ncoef,rank] matrix of xi-coefficients
%   G.cvecs   = {rank,rank}[ncoef,1] cell array of vectors ($c$ in Table 2E)
%   G.dblock  = [ncoef,1] vector of G.smat diagonal indices
%   G.lengths = [order,1] vector of spin dimensions
%   G.invbks  = {rank,rank}[order,1] cell array of 1/<inner products>
%   G.kappas  = {rank,rank}[order,1] cell array of <inner products>
%   G.kvals   = [ncoef,1] vector of order indices
%   G.m1block = [nblock,1] vector of nonzero G.smat row indices
%   G.m2block = [nblock,1] vector of nonzero G.smat column indices
%   G.mranges = [order,2] index of ranges of m-indices
%   G.mvals   = {order,1} cell array of vectors of m-indices
%   G.nblock  = [1] maximal dimension of G.smat
%   G.ncoef   = [1] size of the %G$-blocks (see Table 2D)
%   G.order   = [1] number of columns of $\xi$-values (Table 2C)
%   G.rank    = [1] number of block rows and columns (of $G$-blocks in Table 2D)
%   G.smats   = {rank,rank}[ncoef,ncoef] cell array of sparse matrices ($B+D$ in Table 2F)

% ******************************************************************************
% Copyright (C) 2008 John Arthur Sidles.  This copyright is assigned  
% to the Institute for Soldier Healing of Seattle, Washington, USA.
% ------------------------------------------------------------------------------
% This source file is part of the QSEPACK Template Library (release 1.0d).  
% The QSEPACK Template Library is free software: you can redistribute it 
% and/or modify it under the terms of the GNU General Public License 
% as published by the Free Software Foundation, either version 3 of 
% the License, or (at your option) any later version. 
% 
% The QSEPACK Template Library is distributed in the hope that it will  
% be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
% GNU General Public License for more details. 
% 
% You should have received a copy of the GNU General Public License along  
% with the QSE Template Library. If not, see <http://www.gnu.org/licenses/>.
% 
% The author, John A. Sidles, can be contacted through the Institute for 
% Soldier Healing, 4532 48th Av NE, Seattle, Washington, 98105 USA, or
% by email at "jsidles@soldierhealing.org"
% ------------------------------------------------------------------------------
% The QSEPACK template library is maintained using the open-source
% literate programming tool nuweb (v1.1.1), by Preston Briggs.
% ------------------------------------------------------------------------------
% The QSEPACK template library is distributed as a gzipped tarball file 
% "QSEPACK_v1_0d.tar.gz" that includes the master documentation file 
% "QSEPACK_v1_0d.pdf". 
% ******************************************************************************

% ------------------------------------------------------------------------------
% scrap: store the input as a dense matrix and record its size
if ~iscell(xiIn)
  error('QSEPACK:initialize_metric_QSE:argNotCellArray', '%s', ...
    'the argument xiIn of initialize_metric_QSE(xiIn) must be a cell array.');
end

% OK, store the dimensions
%   G.rank    = [1] number of block rows and columns (of $G$-blocks in Table 2D)
%   G.order   = [1] number of columns of $\xi$-values (Table 2C)
[G.rank,G.order] = size(xiIn);

%   G.lengths = [order,1] vector of spin dimensions
G.lengths  = cellfun('length',xiIn(1,:))'; 

% ... and check for gross errors in input
allRowVectors = true;  % newer versions of MATLAB can use cellfun
for j = 1:G.rank       % to vectorize this test
  for k = 1:G.order
    thisSize = size(xiIn{j,k});
    allRowVectors = allRowVectors & isvector(xiIn(j,k)) & (thisSize(1)==1);
  end
end

if ~allRowVectors
  error('QSEPACK:initialize_metric_QSE:notAllRowVectors', '%s\n%s', ...
    'The cell array argument xiIn of initialize_metric_QSE(xiIn) does ', ...
    'not contain only row vectors (i.e., size(xiIn{j,k}) = [1,n]).');
end

if ~isequal( ...
  cellfun('length',xiIn(:,:)), ...
  repmat(G.lengths',[G.rank,1]) ...
)
  error('QSEPACK:initialize_metric_QSE:argDimensionsInconsistent', '%s\n%s', ...
    'The cell array of row vectors argument xiIn of initialize_metric_QSE(xiIn)', ...
    'has inconsistent dimensions (length(xiIn{j,k}) is not independent of j).');
end

%   G.coefs   = [ncoef,rank] matrix of xi-coefficients
G.coefs = cell2mat(xiIn).'; % flatten and store the data

% ------------------------------------------------------------------------------
% scrap: initialize indexes and size parameters
%   G.ncoef   = [1] size of the %G$-blocks (see Table 2D)
G.ncoef = length(G.coefs(:,1)); 

%   G.mranges = [order,2] index of ranges of m-indices
G.mranges = [ones(size(G.lengths)),G.lengths];
for in = 2:G.order 
  G.mranges(in,:) = G.mranges(in,:) + G.mranges(in-1,2);
end

%   G.mvals   = {order,1} cell array of vectors of m-indices
G.mvals = cell([G.order,1]);
for in = 1:G.order 
  G.mvals{in,1} = (G.mranges(in,1):G.mranges(in,2))';
end

%   G.kvals   = [ncoef,1] vector of order indices
G.kvals  = zeros([sum(G.lengths),1]);
for in = 1:G.order 
  G.kvals(G.mranges(in,1):G.mranges(in,end),1) = in;
end

%   G.m1block = [nblock,1] vector of nonzero G.smat row indices
G.m1block = cell([G.order,1]);
for in = 1:G.order
  G.m1block{in} = reshape( ...
    diag(G.mvals{in})*ones(G.lengths(in)),[],1);
end
G.m1block=cell2mat(G.m1block);

%   G.m2block = [nblock,1] vector of nonzero G.smat column indices
G.m2block = cell([G.order,1]);
for in = 1:G.order
  G.m2block{in} = reshape( ...
    ones(G.lengths(in))*diag(G.mvals{in}),[],1);
end
G.m2block=cell2mat(G.m2block);

%   G.dblock  = [ncoef,1] vector of G.smat diagonal indices
%   G.nblock  = [1] maximal dimension of G.smat
G.dblock = find(G.m2block==G.m1block); 
G.nblock = G.dblock(end); % = sum(G.lengths.*G.lengths)

% ------------------------------------------------------------------------------
% scrap: initialize blockwise kappa values
%   G.invbks  = {rank,rank}[order,1] cell array of 1/<inner products>
%   G.kappas  = {rank,rank}[order,1] cell array of <inner products>
G.kappas = cell([G.rank,G.rank]);
G.invbks = cell([G.rank,G.rank]);
tempbks = complex(zeros([G.order,1]));

for in = 1:G.rank
  for jn = in:G.rank
    tempbks(:) = accumarray( ... % accumarray is build-in
      G.kvals,conj(G.coefs(:,in)).*G.coefs(:,jn));
    G.invbks{in,jn} = 1./tempbks;
    G.kappas{in,jn} = prod(tempbks);
  end
end
clear tempbks;

for in = 2:G.rank
  for jn = 1:(in-1)
    G.invbks{in,jn} = conj(G.invbks{jn,in});
    G.kappas{in,jn} = conj(G.kappas{jn,in});
  end
end

% ------------------------------------------------------------------------------
% scrap: initialize blockwise c vectors
%   G.cvecs   = {rank,rank}[ncoef,1] cell array of vectors ($c$ in Table 2E)
G.cvecs = cell(G.rank);
for i1 = 1:G.rank
  for i2 = 1:G.rank
    G.cvecs{i1,i2} = ... 
      sqrt(G.kappas{i1,i2}).*((G.coefs(:,i2)).* ...
        G.invbks{i1,i2}(G.kvals(:)));
  end
end

% ------------------------------------------------------------------------------
% scrap: initialize a blockwise sparse matrix
%   G.smats   = {rank,rank}[ncoef,ncoef] cell array of sparse matrices ($B+D$ in Table 2F)
G.smats = cell([G.rank,G.rank]);
tempVec = complex(zeros([G.nblock,1]));
for i1 = 1:G.rank
  for i2 = i1:G.rank  % nblock+2*ncoef
    tempVec(:) = (-G.cvecs{i1,i2}(G.m1block)).* ...
      (conj(G.cvecs{i2,i1}(G.m2block)));
    tempVec(G.dblock) = tempVec(G.dblock) + ...
      G.kappas{i1,i2}.*G.invbks{i1,i2}(G.kvals);
    G.smats{i1,i2} = sparse( ...
      G.m1block, ...
      G.m2block, ...
      tempVec,   ...
      G.ncoef,   ...
      G.ncoef,   ...
      G.nblock   ...
    );
  end
end
clear tempVec;

for i1 = 2:G.rank
  for i2 = 1:(i1-1)
    G.smats{i1,i2} = G.smats{i2,i1}';
  end
end